#ifndef __bounce_H_
#define __bounce_H_



namespace math 
{
// numeric integraton. 
template <typename value_type, unsigned dim >
	struct state  {};


struct integration_tag {};
struct euler_intergration : public integration_tag {};
struct runge_kutta_integration : public integration_tag {};


template <typename vector_, typename tag_type = euler_intergration>
	struct integration {};


template <typename vector_ >
	struct integration <vector_, euler_intergration>
{
	typedef vector_ vector_type ;
	typedef typename vector_type::value_type  value_type; 
	
	template <typename unary_function>
	vector_type operator () (const vector_type & xn, const value_type & h, unary_function fun)
	{
		return xn + fun(xn) * h;
	}
};


}


namespace math 
{
// physical modeling
namespace pbm 
{

	struct pbm_system_tag {}; 
	struct bounce_ball : public pbm_system_tag {}; 
	struct other_system : public pbm_system_tag {}; 
	
	template <typename variable_type, typename system_tag>
		struct system_state{} ;

	template <typename variable_>
		struct system_state<variable_, bounce_ball>
	{
		typedef variable_ variable_type ; 
		variable_type velocity, position; 
	};

	template <typename variable_type, typename system_tag>
		struct system_policy {};

	template <typename variable_>
		struct system_policy <variable_, bounce_ball>
	{
		typedef variable_ variable_type ; 
		typedef typename variable_type::value_type value_type; 
		
		static variable_type gravity (const value_type & m) {return variable_type (0.0, 0.0, -1.0) * value_type(g) * value_type(m) ;}
		static variable_type air_force  (const variable_type & velocity)  {return - value_type( d) * velocity ;}
		static variable_type wind_force (const variable_type & wind_dir)  {return  value_type(d) * wind_dir ;}
		static value_type vertical_bounce_velocity (const value_type & vn)  {return - value_type(epsi) * vn ;}
		static value_type horizontal_friction_velocity (const value_type & vt)  {return (1-miu) * vt ;}
		static variable_type acceleration (const variable_type & force)  {return force / value_type(mass); }
		static void decompose_velocity (const variable_type & velocity, const variable_type & normal, variable_type & vn, variable_type & vt)
		{
			vn =  ( velocity * normal ) * normal ;
			vt = velocity - vn ; 
		}
		static variable_type change_collide_velocity (const variable_type & velocity, const variable_type & normal)
		{
			variable_type vn, vt ;
			decompose_velocity (velocity, normal, vn, vt);
			vn = - value_type(epsi) * vn ;
			vt = (1- miu) * vt; 
			return vn + vt; 
		}
	public :
		static constexpr value_type g = 9.8 ; // gravity coefficient ;
		static constexpr value_type d = 6.0 ; // d is air friction coefficient ; 
//		static constexpr value_type time_step = 0.01 ; // h is the timestep ;
		static constexpr value_type epsi = 0.90 ; // bouncing coefficient; 0 \leq epsi \leq 1; 
		static constexpr value_type miu  = 0.05; // miu is horizontal friction coefficient; 
		static constexpr value_type mass = 10; // the mass of the object; 

	};

	template <typename system_tag> struct system {}; 
	template <> struct  system <bounce_ball>
	{
		typedef double value_type; 
		typedef math::vector<value_type, 3> vector_type ;
		typedef vector_type point_type ;
		typedef system_state <math::vector <value_type, 3>, bounce_ball> state_type; 
		typedef system_policy <math::vector <value_type, 3> ,  bounce_ball> policy_type; 

		typedef math::HierarchyMesh<math::mesh_traits<double, 3> >   mesh_type ;
		typedef math::mesh_traits <double, 3>::FaceType              face_type ;
		typedef math::Polygon<value_type, 3>                         polygon_type ;
		typedef math::Sphere<value_type, 3>                          sphere_type ;
		
		void initialize (const state_type & s, value_type ts) 
		{ 
			this->m_state = s;
			time_step = ts ;
		}
		
		vector_type get_acceleration(const vector_type & )
		{
			vector_type force (0.0, 0.0, 0.0); 
			force += policy_type::gravity(value_type (policy_type::mass));
			force += policy_type::air_force (this->m_state.velocity) ;
			//force += policy_type::wind_force (vector_type (1.0, 0.0, 0.0));
			return force / value_type (policy_type::mass); 

		}
		/*
		vector_type get_velocity (const vector_type & ) 
		{
			return vector_type(2.0, 2.0, 2.0) ;
		}
		*/
		state_type moving (value_type remain_time ) 
		{
			// collect all force to check the 
			state_type state  = this->m_state;
			math::integration <vector_type, euler_intergration> obj; 
			value_type h = this->time_step; 

			for (unsigned i = 0; i < 0.01/ h; ++i)
			{
				vector_type old_velocity = state.velocity ;
				state.velocity = obj (old_velocity, h, boost::bind (& system<bounce_ball>::get_acceleration, this, _1)) ;
				//state.position = obj (state.position, h, boost::bind (& system<bounce_ball>::get_velocity, this, _1)) ;
				state.position += old_velocity * h; 
			}
			return state; 
		}

		void face_to_polygon (face_type * face, polygon_type & pol)
		{
			std::vector<point_type *> v_set; 
			face->Face2VertexCoord(v_set) ;
			for (unsigned i = 0; i < v_set.size(); ++i) pol.append_vertex ( * v_set[i]);
		}
		state_type collide (const state_type & old_state, const state_type & new_state, value_type & remain_time) 
		{
			face_type * collide_face = mp_boundary->faces[v_face_id[0]] ;
			vector_type face_normal = collide_face->get_normal (); 
			face_normal.normalize(); 
			polygon_type pol ; 
			face_to_polygon (collide_face, pol); 
			value_type old_dis = std::abs(distance(old_state.position, pol));
			value_type new_dis = std::abs(distance(new_state.position, pol));
			value_type f = old_dis / (old_dis + new_dis);

			state_type state ;
			vector_type collide_velocity = old_state.velocity + f * time_step * get_acceleration(old_state.velocity); 
			state.velocity = policy_type::change_collide_velocity (collide_velocity, face_normal );
			state.position = old_state.position + f * time_step * old_state.velocity; 
			remain_time = time_step * (1-f); 
			return state; 
			
		}
		bool is_collide (const state_type & old_state, const state_type & new_state) 
		{
			v_face_id.clear(); 
			sphere_type s (new_state.position, mp_sphere->get_radius());
			return mp_boundary->intersect(s, v_face_id); 
		}

		void run () 
		{
			static unsigned count = 0; 
			value_type remain_time = time_step; 
			std::cout<<"frame  "<<count++<<std::endl;
			while (remain_time > 0 && m_state.velocity.module2() > 5e-4 )
			{
				state_type new_state = moving(remain_time); 
				if (is_collide (this->m_state, new_state))
				{
					this->m_state = collide (this->m_state, new_state, remain_time) ;
				}
				else
				{
					this->m_state = new_state; 
					remain_time = 0; 
					std::cout<<new_state.velocity<<' '<<new_state.position<<std::endl;
					break; 
				}
			}
		}
		value_type get_time_step() const {return this->time_step ;}
		state_type & get_state () {return this->m_state ;}
		void set_object (mesh_type & mesh, sphere_type & s)
		{
			mp_boundary = & mesh; 
			mp_sphere   = & s; 
		}
	protected:
		
	private:
		state_type m_state; 
		value_type time_step ;
		static constexpr value_type time_max = 10.0; 
		
		mesh_type * mp_boundary; 
		sphere_type * mp_sphere; 

		std::vector<unsigned> v_face_id; 
	};



}
}
#endif
